// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                            avtSPCTHFileFormat.C                           //
// ************************************************************************* //

#include <avtSPCTHFileFormat.h>

#include <string>

#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>
#include <vtkAppendFilter.h>
#include <vtkGeometryFilter.h>
#include <vtkCellData.h>
#include <vtkDataArray.h>

#include <avtDatabaseMetaData.h>
#include <avtMaterial.h>
#include <avtStructuredDomainBoundaries.h>
#include <avtStructuredDomainNesting.h>
#include <avtVariableCache.h>
#include <avtCallback.h>
#include <avtMixedVariable.h>

#include <Expression.h>
#include <InvalidVariableException.h>
#include <InvalidDBTypeException.h>
#include <BadIndexException.h>
#include <DebugStream.h>

using     std::string;


// ****************************************************************************
//  Method: avtSPCTHFileFormat constructor
//
//  Programmer: griffin28 -- generated by xml2avt
//  Creation:   Wed Jun 14 10:05:56 PDT 2017
//
// ****************************************************************************

avtSPCTHFileFormat::avtSPCTHFileFormat(const char * filename) : avtMTMDFileFormat(filename),
                        m_initialized(false), m_fin(NULL), m_spyfile(NULL), m_files(NULL),
                        m_fileCnt(0), m_blocks(0), m_grids(NULL)
{
    m_files = new std::vector<string>();
    m_files->push_back(string(filename));
    
    m_fileCnt = 1;
}

// ****************************************************************************
//  Method: ~avtSPCTHFileFormat destructor
//
//  Programmer: griffin28 -- generated by xml2avt
//  Creation:   Wed Jun 14 10:05:56 PDT 2017
//
// ****************************************************************************

avtSPCTHFileFormat::~avtSPCTHFileFormat()
{
    if(m_initialized)
    {
        for(int i=0; i<m_fileCnt; i++)
        {
            if(m_fin[i] != NULL)
            {
                delete m_fin[i];
            }
            
            if(m_spyfile[i] != NULL)
            {
                delete m_spyfile[i];
            }
        }
    }
    
    if(m_files != NULL)
    {
        m_files->clear();
        delete m_files;
    }
    
    if(m_grids != NULL)
    {
        m_grids->clear();
        delete m_grids;
    }
}

// ****************************************************************************
//  Method: avtSPCTHFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: griffin28 -- generated by xml2avt
//  Creation:   Wed Jun 14 10:05:56 PDT 2017
//
// ****************************************************************************

void
avtSPCTHFileFormat::FreeUpResources(void)
{
    debug4 << "avtSPCTHFileFormat::FreeUpResources" << endl;
}

// ****************************************************************************
//  Method: avtSPCTHFileFormat::Initialize
//
//  Purpose:
//     Read the SPCTH file and check its contents (looks for the magic string
//     "spydata").
//
//
//  Programmer: Kevin Griffin
//  Creation:   Fri June 16 2017
//
//  Modifications:
//
//
// ****************************************************************************

void
avtSPCTHFileFormat::Initialize(void)
{
    if(!m_initialized)
    {
        m_fin = new std::ifstream*[m_fileCnt];
        m_spyfile = new SpyFile*[m_fileCnt];
        
        for(int i=0; i<m_fileCnt; i++)
        {
            m_fin[i] = new ifstream;
            m_fin[i]->open(m_files->at(i).c_str(), ios::in|ios::binary);
            m_spyfile[i] = new SpyFile(m_fin[i]);
            
            char errmsg[500];
            
            // File Header
            if(!m_spyfile[i]->ReadFileHeader())
            {
                snprintf(errmsg, 500, "Error opening SPCTH file %s", GetFilename());
                EXCEPTION1( InvalidDBTypeException, errmsg );
            }
        }
        
        m_initialized = true;
    }
}

// ****************************************************************************
//  Method: avtSPCTHFileFormat::ActivateTimestep
//
//  Purpose:
//      Defines an implementation of ActivateTimestep
//
//  Programmer: Kevin Griffin
//  Creation:   Fri June 16 2017
//
//  Modifications:
//
//
// ****************************************************************************

void
avtSPCTHFileFormat::ActivateTimestep(int ts)
{
    Initialize();
    
    for(int i=0; i<m_fileCnt; i++)
    {
        m_spyfile[i]->ReadTimestepData(ts);
    }
}

// ****************************************************************************
//  Method: avtSPCTHFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: Kevin Griffin -- generated by xml2avt
//  Creation:   Wed Jun 14 10:05:56 PDT 2017
//
// ****************************************************************************

void
avtSPCTHFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md, int ts)
{
    double minExtents[3];
    double maxExtents[3];
    
    for(int i=0; i<m_fileCnt; i++)
    {
        DataBlock *dataBlock = m_spyfile[i]->GetDataBlock(ts);
        
        // Total Blocks
        if(i == 0)
        {
            m_blocks = dataBlock->nBlocks;
        }
        
        // Min/Max Extents
        const double *minE = m_spyfile[i]->GetGmin();
        const double *maxE = m_spyfile[i]->GetGmax();
        
        if(i == 0)
        {
            minExtents[0] = minE[0];
            minExtents[1] = minE[1];
            minExtents[2] = minE[2];
            
            maxExtents[0] = maxE[0];
            maxExtents[1] = maxE[1];
            maxExtents[2] = maxE[2];
        }
        else
        {
            minExtents[0] = minExtents[0] < minE[0] ? minExtents[0] : minE[0];
            minExtents[1] = minExtents[1] < minE[1] ? minExtents[1] : minE[1];
            minExtents[2] = minExtents[2] < minE[2] ? minExtents[2] : minE[2];
            
            maxExtents[0] = maxExtents[0] > maxE[0] ? maxExtents[0] : maxE[0];
            maxExtents[1] = maxExtents[1] > maxE[1] ? maxExtents[1] : maxE[1];
            maxExtents[2] = maxExtents[2] > maxE[2] ? maxExtents[2] : maxE[2];
        }
    }
        
    if(m_blocks > 0)
    {
        //
        //  CODE TO ADD A MESH
        //
        avtMeshMetaData *meshMetaData = new avtMeshMetaData;
        meshMetaData->name = "Mesh";
        meshMetaData->spatialDimension = m_spyfile[0]->GetNumberOfDims();
        meshMetaData->topologicalDimension = m_spyfile[0]->GetNumberOfDims();
        meshMetaData->meshType =  m_fileCnt > 1 ? AVT_UNSTRUCTURED_MESH : AVT_RECTILINEAR_MESH;
        meshMetaData->numBlocks = m_blocks;
        
        meshMetaData->blockOrigin = 0;
        meshMetaData->blockTitle = "patches";
        meshMetaData->blockPieceName = "patch";
        meshMetaData->numGroups = m_spyfile[0]->GetMaxAMRLevels();
        meshMetaData->groupTitle = "levels";
        meshMetaData->groupPieceName = "level";
        
        meshMetaData->containsExteriorBoundaryGhosts = false;
        meshMetaData->hasSpatialExtents = true;
        
        meshMetaData->minSpatialExtents[0] = minExtents[0];
        meshMetaData->minSpatialExtents[1] = minExtents[1];
        meshMetaData->minSpatialExtents[2] = minExtents[2];
        
        meshMetaData->maxSpatialExtents[0] = maxExtents[0];
        meshMetaData->maxSpatialExtents[1] = maxExtents[1];
        meshMetaData->maxSpatialExtents[2] = maxExtents[2];
        
        std::vector<int> groupIds(m_blocks);
        std::vector<std::string> blockPieceNames(m_blocks);
        int levelsOfDetail = 0;
        int idx = 0;
        
        DataBlock *dataBlock = m_spyfile[0]->GetDataBlock(ts);
        int blockCnt = dataBlock->nBlocks;
        
        for(int j=0; j<blockCnt; j++)
        {
            int level = dataBlock->level[j];
            groupIds[idx] = level;
            char temp[128];
            snprintf(temp, 128, "level%d,patch%d", level, idx);
            blockPieceNames[idx] = temp;
            levelsOfDetail = std::max(levelsOfDetail, level);
            
            ++idx;
        }
        
        meshMetaData->blockNames = blockPieceNames;
        meshMetaData->LODs = levelsOfDetail;
        
        md->Add(meshMetaData);
        md->AddGroupInformation(m_spyfile[0]->GetMaxAMRLevels(), m_blocks, groupIds);
        
        //
        // CODE TO ADD A SCALAR VARIABLE
        //
        // string mesh_for_this_var = meshname; // ??? -- could be multiple meshes
        // string varname = ...
        //
        // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
        // avtCentering cent = AVT_NODECENT;
        //
        //
        // Here's the call that tells the meta-data object that we have a var:
        //
        // AddScalarVarToMetaData(md, varname, mesh_for_this_var, cent);
        //
        std::vector<int> matCache;
        const DumpHeader * const dumpHeader = m_spyfile[0]->GetDumpHeader(ts);
        if(dumpHeader != NULL)
        {
            for(int i=0; i<dumpHeader->numVars; i++)
            {
                int fieldId = dumpHeader->vars[i];
                debug5 << "meta::varId = " << fieldId << endl;
                
                // Cell Fields and Derived Fields
                avtScalarMetaData *smd = new avtScalarMetaData;
                if((fieldId > STM_VOID && fieldId < STM_VOLM) || (fieldId >= STM_VM && fieldId <= STM_ZZDEV))
                {
                    smd->meshName = "Mesh";
                    smd->name = m_spyfile[0]->GetFieldName(fieldId);
                    smd->hasUnits = true;
                    smd->units = m_spyfile[0]->GetFieldComment(fieldId);
                    smd->centering = AVT_ZONECENT; //AVT_NODECENT;
                }
                else if((fieldId != STM_VOID) && ((fieldId / STM_VOLM) > 1) )   // Material Fields
                {
                    int matId = m_spyfile[0]->GetMat(fieldId);
                    int matFieldId = fieldId;
                    bool found = false;
                    
                    if(matId >= 0)
                    {
                        matFieldId = fieldId - (matId + 1);
                        
                        for(int i=0; i<matCache.size(); i++)
                        {
                            if(matFieldId == matCache[i])
                            {
                                found = true;
                                break;
                            }
                        }
                    }
                    
                    if(!found)
                    {
                        smd->meshName = "Mesh";
                        smd->name = m_spyfile[0]->GetFieldName(matFieldId);
                        smd->hasUnits = true;
                        smd->units = m_spyfile[0]->GetFieldComment(fieldId);
                        smd->centering = AVT_ZONECENT;
                        matCache.push_back(matFieldId);
                    }
                }
                
                md->Add(smd);
            }
        }
        
        //
        // CODE TO ADD A VECTOR VARIABLE
        //
        // string mesh_for_this_var = meshname; // ??? -- could be multiple meshes
        // string varname = ...
        // int vector_dim = 2;
        //
        // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
        // avtCentering cent = AVT_NODECENT;
        //
        //
        // Here's the call that tells the meta-data object that we have a var:
        //
        // AddVectorVarToMetaData(md, varname, mesh_for_this_var, cent,vector_dim);
        //
        
        //
        // CODE TO ADD A TENSOR VARIABLE
        //
        // string mesh_for_this_var = meshname; // ??? -- could be multiple meshes
        // string varname = ...
        // int tensor_dim = 9;
        //
        // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
        // avtCentering cent = AVT_NODECENT;
        //
        //
        // Here's the call that tells the meta-data object that we have a var:
        //
        // AddTensorVarToMetaData(md, varname, mesh_for_this_var, cent,tensor_dim);
        //
        
        //
        // CODE TO ADD A MATERIAL
        //
        // string mesh_for_mat = meshname; // ??? -- could be multiple meshes
        // string matname = ...
        // int nmats = ...;
        // vector<string> mnames;
        // for (int i = 0 ; i < nmats ; i++)
        // {
        //     char str[32];
        //     sprintf(str, "mat%d", i);
        //     -- or --
        //     strcpy(str, "Aluminum");
        //     mnames.push_back(str);
        // }
        //
        // Here's the call that tells the meta-data object that we have a mat:
        //
        // AddMaterialToMetaData(md, matname, mesh_for_mat, nmats, mnames);
        //
        int numMaterials = m_spyfile[0]->GetNumberOfMats();
        
        if(numMaterials)
        {
            std::vector<string> matNames(numMaterials+1);
            char tmp[32];
            
            for(int i=0; i<numMaterials; i++)
            {
                snprintf(tmp, 32, "mat%d", i+1);
                matNames[i] = tmp;
            }
            
            matNames[numMaterials] = "void";
            
            AddMaterialToMetaData(md, "materials", "Mesh", numMaterials, matNames);
        }
        
        //
        // Here's the way to add expressions:
        //Expression momentum_expr;
        //momentum_expr.SetName("momentum");
        //momentum_expr.SetDefinition("{u, v}");
        //momentum_expr.SetType(Expression::VectorMeshVar);
        //md->AddExpression(&momentum_expr);
        //Expression KineticEnergy_expr;
        //KineticEnergy_expr.SetName("KineticEnergy");
        //KineticEnergy_expr.SetDefinition("0.5*(momentum*momentum)/(rho*rho)");
        //KineticEnergy_expr.SetType(Expression::ScalarMeshVar);
        //md->AddExpression(&KineticEnergy_expr);
        //
    }
}


// ****************************************************************************
//  Method: avtSPCTHFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      ts          The index of the timestate. If GetNTimesteps returned 'N'
//                  times steps, this is guaranteed to be between 0 and N-1.
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: Kevin Griffin -- generated by xml2avt
//  Creation:   Wed Jun 14 10:05:56 PDT 2017
//
// ****************************************************************************

vtkDataSet *
avtSPCTHFileFormat::GetMesh(int ts, int domain, const char *meshname)
{
    if(domain < 0 || domain >= m_blocks)
    {
        EXCEPTION2(BadIndexException, domain, 1);
    }
    
    // Get Mesh dimensions
    bool is3D = m_spyfile[0]->GetNumberOfDims() == 3;
    vtkRectilinearGrid *rg = NULL;
    vtkAppendFilter *appender = NULL;
    
    if(m_grids == NULL)
    {
        m_grids = new std::vector<vtkRectilinearGrid *>();
    }
    
    if(m_fileCnt > 1)
    {
        appender = vtkAppendFilter::New();
    }
    
    for(int fidx=0; fidx<m_fileCnt; fidx++)
    {
        const DataBlock *const dataBlock = m_spyfile[fidx]->GetDataBlock(ts);
        
        if(dataBlock != NULL)
        {
            int dims[3] = {1, 1, 1};
            dims[0] = dataBlock->nx[domain] + 1;
            dims[1] = dataBlock->ny[domain] + 1;
            dims[2] = is3D ? dataBlock->nz[domain] + 1 : 1;
            
            rg = vtkRectilinearGrid::New();
            rg->SetDimensions(dims);
            
            // X Coordinates
            vtkDoubleArray *xcoord = vtkDoubleArray::New();
            xcoord->SetNumberOfTuples(dims[0]);
            xcoord->SetNumberOfComponents(1);
            
            const double *const x = dataBlock->x[domain];
            
            for(int i=0; i<dims[0]; i++)
            {
                xcoord->SetTuple1(i, x[i]);
            }
            
            // Y Coordinates
            vtkDoubleArray *ycoord = vtkDoubleArray::New();
            ycoord->SetNumberOfTuples(dims[1]);
            ycoord->SetNumberOfComponents(1);
            
            const double *const y = dataBlock->y[domain];
            
            for(int i=0; i<dims[1]; i++)
            {
                ycoord->SetTuple1(i, y[i]);
            }
            
            // Z Coordinates
            vtkDoubleArray *zcoord = vtkDoubleArray::New();
            zcoord->SetNumberOfTuples(dims[2]);
            zcoord->SetNumberOfComponents(1);
            
            if(is3D)
            {
                const double *const z = dataBlock->z[domain];
                
                for(int i=0; i<dims[2]; i++)
                {
                    zcoord->SetTuple1(i, z[i]);
                }
            }
            else
            {
                zcoord->SetTuple1(0, 0);
            }
            
            rg->SetXCoordinates(xcoord);
            rg->SetYCoordinates(ycoord);
            rg->SetZCoordinates(zcoord);
            
            xcoord->Delete();
            ycoord->Delete();
            zcoord->Delete();
            
            m_grids->push_back(rg);
            
            if(appender != NULL)
            {
                appender->AddInputData(rg);
            }
        }
    }
    
    if(appender != NULL)
    {
        appender->Update();
        vtkUnstructuredGrid *dset = vtkUnstructuredGrid::New();
        dset->ShallowCopy(appender->GetOutput());
        appender->Delete();
        
        return dset;
    }
    
    return rg;
}

// ****************************************************************************
//  Method: avtSPCTHFileFormat::GetVTKVar
//
//  Purpose:
//      Place the SPY/CTH variable data into a VTK array.
//
//  Arguments:
//      ts          The timestep
//      domain      The current block
//      varname     The SPY/CTH variable name
//
//  Returns:    vtkDataArray
//
//  Programmer: Kevin Griffin
//  Creation:   Wed Jun 14 10:05:56 PDT 2017
//
//  Modifications:
//
// ****************************************************************************

vtkDataArray *
avtSPCTHFileFormat::GetVTKVar(const int ts, const int domain, const char *varname)
{
    if(domain < 0 || domain >= m_blocks)
    {
        EXCEPTION2(BadIndexException, domain, 1);
    }
    
    vtkDoubleArray *scalars = NULL;
    vtkAppendFilter *appender = NULL;
    
    if(m_fileCnt > 1)
    {
        appender = vtkAppendFilter::New();
    }
    
    for(int i=0; i<m_fileCnt; i++)
    {
        DataBlock *const dataBlock = m_spyfile[i]->GetDataBlock(ts);
        int idx = 0;
        
        if(dataBlock != NULL)
        {
            const int fieldId = m_spyfile[i]->GetFieldId(varname);
            debug5 << "fieldId: " << fieldId << endl;
            
            int nx = dataBlock->nx[domain];
            int ny = dataBlock->ny[domain];
            int nz = dataBlock->nz[domain];
            int ncells = nx*ny*nz;
            
            scalars = vtkDoubleArray::New();
            scalars->SetName(varname);
            
            if(ncells > 0)
            {
                scalars->SetNumberOfTuples(ncells);
                
                if(fieldId != -1)
                {
                    for(int k=0; k<nz; k++)
                    {
                        for(int j=0; j<ny; j++)
                        {
                            double *fieldVal = m_spyfile[i]->GetFieldVal(dataBlock, fieldId, k, j, 0, nx, domain);
                            
                            for(int i=0; i<nx; i++)
                            {
                                if(fieldVal != NULL)
                                {
                                    scalars->SetTuple1(idx++, fieldVal[i]);
                                }
                                else
                                {
                                    scalars->SetTuple1(idx++, 0);
                                }
                            }
                        }
                    }
                }
                else
                {
                    for(int i=0; i<ncells; i++)
                    {
                        scalars->SetTuple1(i, 0);
                    }
                }
            }
            else
            {
                scalars->SetNumberOfTuples(1);
                scalars->SetTuple1(1, 0);
            }
            
            if(appender != NULL)
            {
                vtkRectilinearGrid *grid = m_grids->at(i);
                grid->GetCellData()->AddArray(scalars);
                
                appender->AddInputData(grid);
            }
        }
    }
    
    if(appender != NULL)
    {
        appender->Update();
        vtkDoubleArray *data = vtkDoubleArray::New();
        data->DeepCopy(appender->GetOutput()->GetCellData()->GetArray(varname));
        
        // Clean-up
        appender->Delete();
        
        for(int i=0; i<m_fileCnt; i++)
        {
            vtkRectilinearGrid *grid = m_grids->at(i);
            grid->GetCellData()->RemoveArray(varname);
        }
        
        return data;
    }
    
    return scalars;
}

// ****************************************************************************
//  Method: avtSPCTHFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      ts         The index of the timestate. If GetNTimesteps returned 'N'
//                 times steps, this is guaranteed to be between 0 and N-1.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: griffin28 -- generated by xml2avt
//  Creation:   Wed Jun 14 10:05:56 PDT 2017
//
// ****************************************************************************

vtkDataArray *
avtSPCTHFileFormat::GetVar(int ts, int domain, const char *varname)
{
    debug5 << "geting variable " << varname << " for domain: " << domain << " for timestep: " << ts << endl;
    //YOU MUST IMPLEMENT THIS

    //
    // If you have a file format where variables don't apply (for example a
    // strictly polygonal format like the STL (Stereo Lithography) format,
    // then uncomment the code below.
    //
    // EXCEPTION1(InvalidVariableException, varname);
    //

    //
    // If you do have a scalar variable, here is some code that may be helpful.
    //
    // int ntuples = XXX; // this is the number of entries in the variable.
    // vtkFloatArray *rv = vtkFloatArray::New();
    // rv->SetNumberOfTuples(ntuples);
    // for (int i = 0 ; i < ntuples ; i++)
    // {
    //      rv->SetTuple1(i, VAL);  // you must determine value for ith entry.
    // }
    //
    // return rv;
    //
    
    if(domain < 0 || domain >= m_blocks)
    {
        EXCEPTION2(BadIndexException, domain, 1);
    }
    
    const int fieldId = m_spyfile[0]->GetFieldId(varname);
    
    if(fieldId != -1)
    {
        if((fieldId / STM_VOLM) > 1) // isMaterialField(fieldId)
        {            
            int ncells = GetCellCount(ts, domain);
            vtkDoubleArray *scalars = vtkDoubleArray::New();
            
            if(ncells > 0)
            {
                int numMaterials = m_spyfile[0]->GetNumberOfMats();
                std::vector<vtkDoubleArray *> matArrays(numMaterials);
                vtkDataArray *tempArray;
                
                char tmpstr[80];
                for(int m=1; m<=numMaterials; m++)
                {
                    snprintf(tmpstr, 80, "%s+%d", varname, m);
                    tempArray = GetVTKVar(ts, domain, tmpstr);
                    
                    if(tempArray != NULL)
                    {
                        ncells = tempArray->GetSize();
                        matArrays[m-1] = dynamic_cast<vtkDoubleArray *>(tempArray);
                    }
                    else
                    {
                        matArrays[m-1] = NULL;
                    }
                }
                
                scalars->SetNumberOfTuples(ncells);
                // Calculate the cell value
                for(int i=0; i<ncells; i++)
                {
                    double tmpVal = 0;
                    for(int j=0; j<numMaterials; j++)
                    {
                        if(matArrays[j] != NULL)
                        {
                            tmpVal += matArrays[j]->GetValue(i);
                            scalars->SetTuple1(i, tmpVal);
                        }
                    }
                }
                
                // Add Mixed materials
                void_ref_ptr vr = cache->GetVoidRef("materials", AUXILIARY_DATA_MATERIAL, ts, domain);
                avtMaterial *mat = (avtMaterial*) *vr;
                
                if(!mat)
                {
                    DestructorFunction df;
                    void *p = GetAuxiliaryData("materials", ts, domain, AUXILIARY_DATA_MATERIAL, (void*)0, df);
                    mat = (avtMaterial*) p;
                    void_ref_ptr vrtmp = void_ref_ptr(p, df);
                    cache->CacheVoidRef("materials", AUXILIARY_DATA_MATERIAL, ts, domain, vrtmp);
                }
                
                if(mat)
                {
                    std::vector<float> mixvals;
                    for(int z=0; z<mat->GetNZones(); z++)
                    {
                        std::vector<float> vfracs;
                        mat->GetVolFracsForZone(z, vfracs);
                        for(size_t m = 0; m <vfracs.size(); m++)
                        {
                            if(vfracs[m] > 0 && vfracs[m] < 1)
                            {
                                if(m < numMaterials && (matArrays[m] != NULL))
                                {
                                    mixvals.push_back(matArrays[m]->GetTuple1(z));
                                }
                                else
                                {
                                    mixvals.push_back(0.0f);
                                }
                            }
                        }
                    }
                    
                    if ((int) mixvals.size() != mat->GetMixlen())
                    {
                        char errmsg[256];
                        snprintf(errmsg, sizeof(errmsg), "Mixed variable size, %d, doesn't agree with material "
                                 "mixlen, %d, for variable \"%s\"", (int) mixvals.size(), mat->GetMixlen(), varname);
                        //                    avtCallback::IssueWarning(errmsg);
                        debug5 << errmsg << endl;
                    }
                    else
                    {
                        avtMixedVariable *mv = new avtMixedVariable(&mixvals[0], (int) mixvals.size(), varname);
                        vr = void_ref_ptr(mv, avtMixedVariable::Destruct);
                        cache->CacheVoidRef(varname, AUXILIARY_DATA_MIXED_VARIABLE, ts, domain, vr);
                    }
                }
                
                // Clean-up
                for(int i=0; i < numMaterials; i++)
                {
                    if(matArrays[i] != NULL)
                    {
                        matArrays[i]->Delete();
                    }
                }
            }
            else
            {
                scalars->SetNumberOfTuples(1);
                scalars->SetTuple1(0, 0);
            }
            
            return scalars;
        }
        else    // Cell Field
        {
            return GetVTKVar(ts, domain, varname);
        }
    }
    else
    {
        EXCEPTION1(InvalidVariableException, varname);
    }
    
    return NULL;
}


// ****************************************************************************
//  Method: avtSPCTHFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      ts         The index of the timestate. If GetNTimesteps returned 'N'
//                 times steps, this is guaranteed to be between 0 and N-1.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: griffin28 -- generated by xml2avt
//  Creation:   Wed Jun 14 10:05:56 PDT 2017
//
// ****************************************************************************

vtkDataArray *
avtSPCTHFileFormat::GetVectorVar(int ts, int domain, const char *varname)
{
    // YOU MUST IMPLEMENT THIS
    //
    // If you have a file format where variables don't apply (for example a
    // strictly polygonal format like the STL (Stereo Lithography) format,
    // then uncomment the code below.
    //
    // EXCEPTION1(InvalidVariableException, varname);
    //

    //
    // If you do have a vector variable, here is some code that may be helpful.
    //
    // int ncomps = YYY;  // This is the rank of the vector - typically 2 or 3.
    // int ntuples = XXX; // this is the number of entries in the variable.
    // vtkFloatArray *rv = vtkFloatArray::New();
    // int ucomps = (ncomps == 2 ? 3 : ncomps);
    // rv->SetNumberOfComponents(ucomps);
    // rv->SetNumberOfTuples(ntuples);
    // float *one_entry = new float[ucomps];
    // for (int i = 0 ; i < ntuples ; i++)
    // {
    //      int j;
    //      for (j = 0 ; j < ncomps ; j++)
    //           one_entry[j] = ...
    //      for (j = ncomps ; j < ucomps ; j++)
    //           one_entry[j] = 0.;
    //      rv->SetTuple(i, one_entry); 
    // }
    //
    // delete [] one_entry;
    // return rv;
    //
    return 0;
}

// ****************************************************************************
//  Method: avtSPCTHFileFormat::GetAuxiliaryData
//
//  Purpose:
//      Gets the auxiliary data specified.
//
//  Arguments:
//      var         The variable of interest.
//      ts          The timestep
//      domain      The domain of interest.
//      type        The type of auxiliary data.
//      args        The arguments for that type.
//
//  Returns:    The auxiliary data.  Throws an exception if this is not a
//              supported data type.
//
//  Programmer: Kevin Griffin
//  Creation:    Wed Jun 14 10:05:56 PDT 2017
//
// ****************************************************************************

void *
avtSPCTHFileFormat::GetAuxiliaryData(const char *var, int ts, int domain, const char *type,
                                     void *args, DestructorFunction &df)
{
    void *rv = NULL;
    
    if(strcmp(type, AUXILIARY_DATA_MATERIAL) == 0)
    {
        rv = GetMaterial(var, ts, domain);
        df = avtMaterial::Destruct;
    }
    
    return rv;
}

// ****************************************************************************
//  Method: avtSPCTHFileFormat::GetMaterial
//
//  Purpose:
//      Gets the material auxiliary data.
//
//  Arguments:
//      var         The variable of interest.
//      ts          The timestep
//      domain      The domain of interest.
//
//  Returns:    The material data
//
//  Programmer: Kevin Griffin
//  Creation:    Wed Jun 14 10:05:56 PDT 2017
//
// ****************************************************************************

void *
avtSPCTHFileFormat::GetMaterial(const char *var, const int ts, const int domain)
{
    debug5 << "Reading in timestep " << ts << " domain " << domain << ", material " << var << endl;
    
    if(domain < 0 || domain >= m_blocks)
    {
        EXCEPTION2(BadIndexException, domain, 1);
    }
    
    int numMaterials = m_spyfile[0]->GetNumberOfMats();
    std::vector<string> matNames(numMaterials+1);
    char tmp[32];
    
    for(int i=0; i<numMaterials; i++)
    {
        snprintf(tmp, 32, "mat%d", i+1);
        matNames[i] = tmp;
    }
    
    matNames[numMaterials] = "void";
    int ncells = GetCellCount(ts, domain);
    
    // Mixed Materials
    std::vector<int>    mix_zone;
    std::vector<int>    mix_mat;
    std::vector<float>  mix_vf;
    std::vector<int>    mix_next;
    
    // Get the material fractions
    std::vector<vtkDoubleArray *> matArrays(numMaterials+1);
    vtkDataArray *tempArray;
    
    // Get Material Fractions
    char tmpstr[80];
    for(int i=1; i<=numMaterials; i++)
    {
        snprintf(tmpstr, 80, "VOLM+%d", i);
        tempArray = GetVTKVar(ts, domain, tmpstr);
        
        if(tempArray != NULL)
        {
            ncells = tempArray->GetSize();
            matArrays[i-1] = dynamic_cast<vtkDoubleArray *>(tempArray);
        }
        else
        {
            matArrays[i-1] = NULL;
        }
    }
    
    std::vector<int>    matlist(ncells);
    
    // Get Void Fractions
    tempArray = GetVTKVar(ts, domain, m_spyfile[0]->GetFieldName(STM_VOID));
    
    if(tempArray != NULL)
    {
        matArrays[numMaterials] = dynamic_cast<vtkDoubleArray *>(tempArray);
    }
    else
    {
        matArrays[numMaterials] = NULL;
    }
    
    // Build the mixed material data structures
    for(int i=0; i<ncells; i++)
    {
        // Pure materials
        int numCellMats = 0;
        int lastCellMat = -1;
        
        for(int j=0; j<=numMaterials; j++)
        {
            if(matArrays[j] != NULL)
            {
                if(matArrays[j]->GetValue(i) > 0)
                {
                    numCellMats++;
                    lastCellMat = j;
                }
            }
        }
        
        if(numCellMats == 1)
        {
            matlist[i] = lastCellMat;
            continue;
        }
        
        if(numCellMats == 0)
        {
            matlist[i] = numMaterials;
            continue;
        }
        
        // Mixed Materials
        bool addRemainder = false;
        double vfRemainder = 1;
        
        for(int j=0; j<=numMaterials; j++)
        {
            vfRemainder -= matArrays[j]->GetValue(i);
        }

        if(std::abs(vfRemainder) > 1.0e-5)
        {
            addRemainder = true;
        }
        matlist[i] = -1 * (1 + mix_zone.size());
        for(int j=0; j<=numMaterials; j++)
        {
            double matValue = matArrays[j]->GetValue(i);
            if(addRemainder && (j == numMaterials))
            {
                matValue += vfRemainder;
            }
            
            debug5 << "Material = " << j << "Value = " << matValue << endl;
            
            if(matValue <= 0)
            {
                continue;
            }
            
            mix_zone.push_back(i);
            mix_mat.push_back(j);
            mix_vf.push_back(matValue);
            mix_next.push_back(mix_zone.size()+1);
        }
        
        // Last entry is zero
        if(mix_next.size() > 0)
        {
            mix_next[mix_next.size()-1] = 0;
        }
    }
    
    // Clean-up
    for(int i=0; i <= numMaterials; i++)
    {
        if(matArrays[i] != NULL)
        {
            matArrays[i]->Delete();
        }
    }
    
    int mixed_size = mix_zone.size();
    int *ml = NULL, *mixm = NULL, *mixn = NULL, *mixz = NULL;
    float *mixv = NULL;
    
    if(matlist.size() > 0)
    {
        ml = &(matlist[0]);
    }
    
    if(mix_zone.size() > 0)
    {
        mixz = &(mix_zone[0]);
    }
    
    if(mix_mat.size() > 0)
    {
        mixm = &(mix_mat[0]);
    }
    
    if(mix_vf.size() > 0)
    {
        mixv = &(mix_vf[0]);
    }
    
    if(mix_next.size() > 0)
    {
        mixn = &(mix_next[0]);
    }
    
    avtMaterial *mat = new avtMaterial(numMaterials, matNames, ncells, ml, mixed_size, mixm, mixn, mixz, mixv);
    
    return (void *)mat;
}

// ****************************************************************************
//  Method: avtSPCTHFileFormat::GetNTimesteps
//
//  Purpose:
//      Get the number of timesteps
//
//  Returns:    The number of timesteps
//
//  Programmer: Kevin Griffin
//  Creation:    Wed Jun 14 10:05:56 PDT 2017
//
// ****************************************************************************

int
avtSPCTHFileFormat::GetNTimesteps()
{
    Initialize();
    return m_spyfile[0]->GetNumberOfDumps();
}

// ****************************************************************************
//  Method: avtSPCTHFileFormat::GetCycles
//
//  Purpose:
//      Get the cycles
//
//
//  Programmer: Kevin Griffin
//  Creation:    Wed Jun 14 10:05:56 PDT 2017
//
// ****************************************************************************

void
avtSPCTHFileFormat::GetCycles(std::vector<int> &cycles)
{
    Initialize();
    const std::vector<Dump *> * const dumps = m_spyfile[0]->GetDumps();
    
    for(int i=0; i<dumps->size(); i++)
    {
        cycles.push_back(dumps->at(i)->cycle);
    }
}

// ****************************************************************************
//  Method: avtSPCTHFileFormat::GetTimes
//
//  Purpose:
//      Get the times
//
//
//  Programmer: Kevin Griffin
//  Creation:    Wed Jun 14 10:05:56 PDT 2017
//
// ****************************************************************************

void
avtSPCTHFileFormat::GetTimes(std::vector<double> &times)
{
    Initialize();
    const std::vector<Dump *> * const dumps = m_spyfile[0]->GetDumps();
    
    for(int i=0; i<dumps->size(); i++)
    {
        times.push_back(dumps->at(i)->dTime);
    }
}

// ****************************************************************************
//  Method: avtSPCTHFileFormat::SetFileList
//
//  Purpose:
//      Stores the list of files contained in the .visit file.
//
//  Arguments:
//      list         list of files from the .visit file.
//      nlist        the number of files in the list
//
//
//  Programmer: Kevin Griffin
//  Creation:   Wed Dec 11 10:05:56 PDT 2017
//
// ****************************************************************************

void
avtSPCTHFileFormat::SetFileList(const char *const *list, const int nList)
{
    string str = string("*.visit");
    
    delete filename;
    filename = new char[strlen(str.c_str())+1];
    strcpy(filename, str.c_str());
    
    if(m_files == NULL)
    {
        m_files = new std::vector<string>();
    }
    else
    {
        m_files->clear();
    }
    
    for(int i=0; i<nList; i++)
    {
        m_files->push_back(string(list[i]));
    }
    
    m_fileCnt = nList;
}

// ****************************************************************************
//  Method: avtSPCTHFileFormat::GetCellCount
//
//  Purpose:
//      Gets the maximum number of cells accross all files for the current
//      timestep and domain.
//
//  Arguments:
//      ts          the current timestep
//      domain      the current block
//
//  Returns:        the number of cells
//
//  Programmer:     Kevin Griffin
//  Creation:    Wed Dec 11 10:05:56 PDT 2017
//
// ****************************************************************************

int
avtSPCTHFileFormat::GetCellCount(const int ts, const int domain)
{
    int xdim = 0;
    int ydim = 0;
    int zdim = 0;
    
    for(int idx=0; idx<m_fileCnt; idx++)
    {
        const DataBlock *const dataBlock = m_spyfile[idx]->GetDataBlock(ts);
        
        if(dataBlock != NULL)
        {
            xdim += dataBlock->nx[domain];
            ydim += dataBlock->ny[domain];
            zdim += dataBlock->nz[domain];
        }
    }
    
    debug5 << "nx: " << xdim << endl;
    debug5 << "ny: " << ydim << endl;
    debug5 << "nz: " << zdim << endl;
    
    int ncells = m_spyfile[0]->GetNumberOfDims() == 2 ? (xdim * ydim) : (xdim * ydim * zdim);
    return ncells;
}


